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Abstract 

The highly convergent iterative methods developed by 
Trujillo Bueno and Fabiani Bendicho (1995) for radia- 
tive transfer (RT) applications are generalized to spher- 
ical symmetry with velocity fields. These RT methods are 
based on Jacobi, Gauss-Seidel (GS), and SOR iteration 
and they form the basis of a new NLTE multilevel trans- 
fer code for atomic and molecular lines. The benchmark 
tests carried out so far are presented and discussed. The 
main aim is to develop a number of powerful RT tools for 
the theoretical interpretation of molecular spectra. 
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1. Introduction 

Atomic line emission has been extensively used for trac- 
ing the physical conditions in many astrophysical plasmas. 
The relatively high energy difference between the atomic 
energy levels makes this diagnostic tool a suitable one for 
tracing the physical conditions in warm and hot media. 
But cool plamas are not well traced by atomic lines be- 
cause the thermal energy is not high enough to populate 
the upper levels of the transitions. Fortunately, FIRST 
will allow us to study molecular line emission in greater 
detail. Molecules have very rich spectra, arising from tran- 
sitions between the electronic, vibrational and rotational 
levels. Their spectra cover the spectral range from radio to 
optical wavelengths, depending on the type of transition. 
One has to include very complicated molecular systems to 
be able to model the observations and an extensive for- 
ward RT modeling effort is frequently needed. With this 
motivation in mind, we are developing a radiative trans- 
fer code based on the fast iterative methods developed 



by Trujillo Bueno & Fabiani Bendicho (1995) for Carte- 
sian coordinates. Our first step was to generalize these 
methods to spherical symmetry with velocity fields. This 
new RT tool will help us in the interpretation of different 
kinds of observations, including ro-vibrational bands in 
circumstellar envelopes of C-rich or 0-rich evolved stars, 
rotational lines in molecular clouds, molecular emission 
from the Sun, maser emission and many others. Another 



interesting application would be the modeling of maser 
polarization. 

2. The state of the art and our approach 

The radiative transfer problem requires the self-consistent 
solution of the rate equations for the populations of the 
molecular levels and the radiative transfer equation. This 
set of equations describes a nonlinear and nonlocal prob- 
lem. Iterative methods are therefore needed. Since Bernes 
1(1979) , most of the radiative transfer tools used in molec- 
ular radiative transfer have been based on Monte Carlo 
techniques for the solution of the radiative transfer equa- 
tion and on the A-iteration for the iterative solution of 
the non-LTE problem. The Monte Carlo technique is very 
powerful for its ability to cope with complicated geome- 
tries, but it has a major drawback: its instrinsic random 
noise. There have been many efforts to reduce the noise. 



but it is always present (see, for example, Bernes 1979). 
On the other hand, the A-iteration scheme is very easy to 
implement, but is only useful for optically thin problems. 
Recently, Accelerated Lambda Iteration (ALI) methods 
have been applied to molecular radiative transfer. ALI is 
a method that requires not much more computational time 
per iteration than the A-iteration, but it has a better con- 
vergence behavior for optically thick problems. It has been 
extensively used in stellar and solar astrophysics and has 
recently been combined with Monte Carlo techniques for 
molecular RT ( Hogerheijde fc van der Tak 2000 ). 



Our approach is based on the iterative methods devel- 
oped by Trujillo Bueno fc Fabiani Bendicho (1995) , which 
themselves are based on the Gauss-Seidel scheme and ap- 
plied to Cartesian coordinates in ID, 2D and 3D, either 
with or without polarization. They allow the solution of 
the radiative transfer problem with the same computa- 
tional time per iteration as ALI, but with an order-of- 
magnitude improvement in the number of iterations. The 
short-characteristics technique ( Kunasz fc Auer 198^ ) has 



been chosen as the formal solver. This has facilitated the 
implementation of the scheme in a very efficient way (see 
Trujillo Bueno fc Fabiani Bendicho 1995| ). We have gener- 
alized these methods to spherical symmetry with velocity 
fields. The problem is still one dimensional because the 
physical variables have only radial dependence, but an- 
gular information has to be achieved more precisely than 
for a plane-parallel atmosphere in order to take account 
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of curvature effects. This angular information is obtained 
by means of solving the RT equation through the im- 
pact parameters, as is usually done. The velocity fields 
are treated in the observer's frame. In Fig. (|l|) we show 
a schematic representation of the difference between the 
ALI-based and the GS-based iterations. With ALI, when 
the radiation field is obtained at all the points of the 
atmosphere, the statistical equilibrium equations for the 
molecular population can be solved. On the other hand, 
the main idea behind the GS-based methods is the fact 
that when the radiation field is known at one point in 
the atmosphere, one can write the statistical equilibrium 
equations for this point and do the level population correc- 
tion. When solving the RT equation to get the radiation 
field at the next point, we have to take into account that 
the population in the previous point has been improved. 
This scheme, coded in an efficient way with the aid of 
the short-characteristics formal solver, can lead to a high 
improvement in the total number of iterations, while the 
time per iteration is virtually the same as with ALI. 



temperature for the transitions J = 1 ^ and ,7 = 2-^1 
through the cloud. We show the results obtained with our 
code and that obtained by Bernes using his Monte Carlo 
technique. The intrinsic noise of the Monte Carlo scheme 
can clearly be appreciated in the figures. 
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Figure 2. Excitation temperature for the J = 1 — > and 
J — 2 1 rotational transitions of CO in the Bernes' 
cloud. Comparison between the our results and those of 
Bernes (1979) are plotted. 
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Figure 1. Scheme showing the differences between the ALI 
iterative method and the GS-based method. The figure is 
explained in the text. 



3. Illustrative examples 

In order to verify that our code is giving reliable results, 
we have chosen several benchmark problems whose results 
have already been published. 

3.1. BeRNES'S CO CLOUD 



Let us begin with the CO cloud model used by Bernes 
(1979) to introduce his Monte Carlo code. The problem 



consists in a constant-density (riHa — 2 x 10 cm^'^), 
constant-temperature (T = 20 K), 1 pc radius infalling 
cloud with a maximum velocity of 1 km s^^ at the external 
parts and sampled at 40 radial shells. The CO abundance 
is 5 X lO^'^ and the cloud is illuminated by the cosmic mi- 
crowave background radiation (CMBR) at a temperature 
of 2.7 K. The CO molecule with the first six rotational lev- 
els is used, taking into account that the same coUisional 
rates used by Bernes in his calculation have to be used 
to get similar results. In Fig. (0) we show the excitation 



3.2. Leiden benchmark test 

A number of useful test cases became available after the 
1999 workshop on Radiative Transfer in Molecular Lines 
at the Lorentz Center of Leiden University]^. These are 
intended for the testing of newly developed molecular RT 
codes against already existing ones. Although every test 
problem has been solved with our multilevel NLTE code 
and good agreement obtained, we show only some of the 
results. The model describes a collapsing cloud similar to 
that described by Shu (1977), where the first 21 rotational 
levels of HCO+ (from J = to J = 20) are taken into 
account in the non-LTE calculation. The molecular abun- 
dance is [HCO^] = 10~^, so lines are only slightly opti- 
cally thick (r < 10). The cloud is sampled logarithmically 
at 50 depth points and is externally illuminated by the 
CMBR at 2.728 K. Results for J = and J = 1 are given 
if Fig. (|3^) for the different codes used in the test and in 
Fig. (||b) corresponding to our code. In these plots we rep- 
resent the fractional population of each level, which can 
be written as / = nievei/^totai- 

We see that the results agree. Although it cannot be 
seen in our plots, most of the HCO"'' in the inner parts of 
the cloud is in the lowest four rotational levels, because 
the kinetic temperature is relatively high {T < 20 K) and 
there is energy in the medium to populate the higher lev- 
els. On the other hand, in the external zones of the cloud, 
almost 60% of the IICO+ is at the ground level (J = 0), 
and the remaining 40% is in the J — 1 level due to the 
lower kinetic temperature, which is not able to populate 
the higher levels efficiently. As can be seen in Fig. (||a), 
there is one curve which is different from the others. This 
is caused by not having included the CMB radiation as the 

^ http:/ /www. strw.leidenuniv.nl/~radtrans 
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Figure 3. Fractional population of the first two rotational 
levels of HCO^ . Panel a) represents the results obtained 
with seven different codes, while panels b) and c) represent 
the results obtained with our code by including (b) and 
excluding (c) the CMBR. 



Figure 4.. Excitation temperature for the rotational tran- 
sitions between the three lowest levels of HCO^ . Panel a) 
represents the results obtained with seven different codes, 
while panels (b) and (c) represent the results obtained with 
our code, either including the CMBR (b) or not (c). 



outer boundary condition and assuming that the cloud is 
not externaUy illuminated. This turns out to be an extra 
test for our code, and the results for this particular situa- 
tion are shown in Fig. (^). Agreement is also obtained for 
the remaining levels, and one can see that there is no sig- 
nificant diff'erence between the results in the inner parts of 
the cloud. However, a totally diff'erent result is obtained in 
the external zones, where the ground level is the only one 
populated with ^--^90% of the total abundance. Although 
the kinetic temperature at the outer envelopes of the cloud 
is still able to populate higher levels, non-LTE effects pro- 
duce this underpopulation of the higher levels. Excitation 
temperature is also plotted in Fig. (^, for the two transi- 
tions J = 1 — !• and the J = 2 ^ 1. Also in Figs. (||b) and 
(^) the results obtained with our code are also plotted, 
either including the CMBR as a boundary condition or 
not, respectively. The results are also comparable to those 
obtained by different codes in both cases. There is a little 
more dispersion in this result than in that for fractional 
population, but this could be due to the fact that the 
majority of the codes are based on Monte Carlo schemes, 
which, although they have variance reduction techniques, 
have an intrinsic random noise that could produce these 
effects. 



3.3. CO IN THE Sun 

We have solved the non-LTE problem for the A?; = 1 ro- 
vibrational band of CO at 4.7 /xm. The vibrational con- 
stant for CO is cjo = 2143 cm^^ and the rotational con- 
stant for the vibrational ground state is Bo — 1.923 cm~^. 
Since ujq ^ By, it follows that successive rotational lev- 
els within one vibrational state are much closer in en- 
ergy than similar rotational levels in successive vibrational 
states. It can be shown that spontaneous radiative decay 
rates for pure rotational transitions are much lower than 



coUisional rates (see, for example, Thompson 1973), so 
one may assume without many problems that the popula- 
tions of the rotational levels within a vibrational state are 
given by Boltzmann statistics. This assumption greatly 
simplifies the problem as shown by Uitenbrock (2000")| for 
the same CO problem, because the number of unknowns 
is reduced from the total number of levels (the number 
of vibrational levels x number of rotational levels within 
each vibrational state) to the total number of vibrational 
levels. However, we have solved the whole problem with- 
out making this assumption and including the first five 
vibrational levels (from v ^ to v ~ A) and 21 rota- 
tional levels within each vibrational one (from J = to 
J = 20). A quiet-Sun model atmosphere has been cho- 
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Figure 5. Departure coefficients ( right axis ) and CO abun- 
dance (left axis) in a non-LTE calculation for the Av = 1 
hand at ^.1 [im in a quiet-Sun model. The first five vi- 
brational levels (v — to v = A) with 21 rotational levels 
(J = to J = 20) within each one are included. Note that 
LTE is obtained for CO in the line-formation region. 



sen ( Vernazza et al. 1981 ) and the molecular abundance 
has been calculated in this model assuming chemical equi- 
librium. As shown in Fig. (H), the CO abundance peaks 
at ^ 300 km above the bottom of the photosphere. This 
figure also shows the departure coefficients for all the vi- 
brational levels included. These departure coefficients are 
calculated as usual, but taking into account the total pop- 
ulation of each vibrational level, which can be obtained 
by summing over the rotational levels inside each vibra- 
tional one: h — X]j "nlte(w, J)/ X]j "lteI'^, >/)■ We see 
that LTE is obtained in the line-formation region below 
~ 800 km, which is the zone where most of the CO is 
formed. This partially confirms the results that can be 
obtained by comparison of the radiative and coUisional 
transition rates. 



the possible important levels is crucial for the interpreta- 
tion of observations. On the other hand, the advantage of 
getting the solution of the non-LTE problem in only a few 
iterations leads to great advantages. This makes it possi- 
ble to improve the adjusting of all the physical parameters 
in the model one is using to interpret the observations, be- 
cause much more extensive forward modeling is now pos- 
sible. Finally, the fast solution of the radiative transfer 
problem allows us to introduce the transfer of polarized 
radiation wit h the aid of the den sity matrix theory (see 
the review by Trujillo Bueno 2001 ) . This could lead to the 
self-consistent solution of the maser polarization problem, 
which could make it possible to explain radio observations 
of masers such as SiO. 
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4. Conclusions 

We have generalized very efficient iterative methods for 
the solution of the molecular radiative transfer problem 
to spherical geometry with velocity fields. The problem 
is still one-dimensional, but more angular information is 
required in comparison to the plane-parallel case, so the 
total computation time is larger. Velocity fields are treated 
in the observer's frame, so velocity fields have to be lim- 
ited to several times the thermal velocity in the medium 
if we want to have a tractable frequency quadrature. This 
limitation is only a computational problem and not a true 
limitation of the method. Such a fast solution of the non- 
LTE problem allows the solution of more complicated sit- 
uations, where larger molecular models can be used. It is 
known that there can be many difi'erent pumping processes 
in molecular radiative transfer — very important for the in- 
terpretation of masers — and a correct model including all 



